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The topological insulator is an electronic phase stabilized by spin-orbit coupling that supports 
propagating edge states and is not adiabatically connected to the ordinary insulator. In several 
ways it is a spin-orbit-induced analogue in time-reversal-invariant systems of the integer quantum 
Hall effect (IQHE) . This paper studies the topological insulator phase in disordered two-dimensional 
systems, using a model graphene Hamiltonian introduced by Kane and Mele as an example. The 
nonperturbative definition of a topological insulator given here is distinct from previous efforts in 
that it involves boundary phase twists that couple only to charge, does not refer to edge states, 
and can be measured by pumping cycles of ordinary charge. In this definition, the phase of a 
Slater determinant of electronic states is determined by a Chern parity analogous to Chern number 
in the IQHE case. Numerically we find, in agreement with recent network model studies, that 
the direct transition between ordinary and topological insulators that occurs in band structures is 
a consequence of the perfect crystalline lattice. Generically these two phases are separated by a 
metallic phase, which is allowed in two dimensions when spin-orbit coupling is present. The same 
approach can be used to study three-dimensional topological insulators. 
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I. INTRODUCTION 

Considerable theoretical and experimental effort has 
bee n devo ted to the quest for an intrinsic spin Hall ef- 
feclP'SlSIll that would allow generation of spin currents by 
an applied electric field. Interesting mechanisms for such 
spin current generation make use of spin-orbit coupling, 
which breaks the SU (2) spin symmetry of free electrons 
but not time-reversal symmetry. A dissipationless type of 
intrinsic spin Hall effect was predicted--^ to arise in mate- 
rials that have an electronic energy gap. This "quantum 
spin Hall effect" (QSHE) in certain materials with time- 
reversal symmetry has a subtle relationship to the integer 
quantum Hall effect, in which time-reversal symmetry is 
explicitly broken by a magnetic field. 

In a system with unbroken time-reversal symmetry, a 
dissipationless charge current is forbidden, but a dissipa- 
tionless transverse spin current is allowed, of the form 

j; = ae,,kEk. (1) 

The current on the left is a spin current and e is the 
fully antisymmetric tensor. Note that a spin current re- 
quires two indices, one for the direction of the current and 
one for the direction of angular momentum that is trans- 
ported. The constant of proportionality a depends on the 
specific mechanism: for example, the (dissipative) extrin- 
sic D 'yakonov-Perel mechanisnl^ predicts a small a that 
depends on impurity concentration. The QSHE builds 
on the construction by HaldancP of a lattice "Chern in- 
sulator" model, with broken time-reversal symmetry but 
without net magnetic flux, that shows a v ^ \ IQHE. 
The simplest example of a QSHE is obtained by taking 
two copies of Haldane's model, one for spin-up electrons 
along some axis and one for spin-down. Time-reversal 
symmetry can be maintained if the effective IQHE mag- 



netic fields are opposite for the two spin components. 
Then an applied electric field generates a transverse cur- 
rent in one direction for spin-up electrons, and in the op- 
posite direction for spin-down electrons. There is no net 
charge current, consistent with time-reversal symmetry, 
but there is a net spin current. However, models like this 
in which one component of spin is perfectly conserved are 
both unphysical, since realistic spin-orbit coupling does 
not conserve any component, and not very novel, since 
for each spin component the physics is exactly the same 
as Haldane's model and the spin components do not mix. 

More subtle physics emerges when one asks how 
the QSHE appears in more realistic band structures. 
Remarkably, band insulators of noninteracting two- 
dimensional electrons with spin-orbit coupling divide into 
two classes: the "ordinary insulator" , which in general 
has no propagating edge modes and no spin Hall effect, 
and the "topological insulator", which has stable prop- 
agating edge modes and a generic spin Hall effect, al- 
though the amount of spin transported (the coefficient 
a in ([l])) is nonuniversal. These phases are associated 
with a Z2-valued topological invariant (an "oddness" 
or "evenness" , in the language of parity* ) in the same 
way that IQHE phases are associated with an integer- 
valued topological invariant. For explicitness, consider 
the model of graphene introduced by Kane and Mele.'^ 
This is a tight-binding model for independent electrons 
on the honeycomb lattice (Fig.jl]). The spin-independent 
part of the Hamiltonian consists of a nearest-neighbor 
hopping, which alone would give a semimetallic spectrum 



* Throughout "parity" is used to denote oddness or evenness, 
rather than a spatial inversion eigenvalue. 
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with Dirac nodes at certain points in the 2D Brillouin 
zone, plus a staggered sublattice potential whose effect is 
to introduce a gap: 



(2) 



Here {ij) denotes nearest-neighbor pairs of sites, cr is a 
spin index, alternates sign between sublattices of the 
honeycomb, and t and Xy are parameters. 

The insulator created by increasing is an unre- 
markable band insulator. However, the symmetries of 
graphene also permit an "intrinsic" spin-orbit coupling 
of the form 

Hso = iXsO ''«j4o-i<<T2Cjt2- (3) 

((ij>)CTi(T2 

Here = {2/\/3)di x d2 = ±1, where i and j are next- 
nearest-neighbors and di and ^2 are unit vectors along 
the two bonds that connect i to j. The Hamiltonian 
Hq + Hso conserves , the distinguished component of 
electron spin, and reduces for fixed spin (up or down) to 
Haldane's model.'^ 




FIG. 1: The honeycomb lattice on which the tight-binding 
Hamiltonian resides. For the two sites depicted, the factor 



of equation Q 



= —1. The phases 



describe 



twisted boundary conditions, introduced in equation ( 10 1 



The "topological insulator" phase created when 
|Aso| ^ |A«| is quite different from the ordinary insu- 
lator that appears when |A„| ^ \Xso\ (here we assume 
that there is an energy gap between the lower and up- 
per band pairs in which the Fermi level lies). The for- 
mer has counterpropagating edge modes and shows the 
QSHE, while the latter does not.'^Does this phase exist 
for more realistic spin-orbit coupling? The spin compo- 
nent Sz is no longer a good quantum number when the 
Rashba spin-orbit coupling is added: 



iX 



(72 ^ ^ij 



(4) 



with dij the vector from i ^ j and dij the correspond- 
ing unit vector. (Note that Rashba spin-orbit coupling 
is not intrinsic to graphene but generated by inversion- 
symmetry breaking in the out-of-plane direction.) The 
topological insulator survives but is strongly modified in 
the presence of this term. For a general 2D band struc- 
ture with conserved, there are many phases labeled 
by an integer n, as in the IQHE: if spin-up electrons are 
in the ly ~ n state, then spin-down electrons must be in 
the v = —n state by time-reversal symmetry, where the 
sign indicates that the direction of the effective magnetic 
field is reversed. Once is not conserved, there are only 
two insulating phases, the "ordinary" and "topological" 
insulators. A heuristic definition of the topological insu- 
lator, without reference to any particular spin component 
or the spin Hall effect, is as a band insulator that is re- 
quired to have gapless propagating edge modes at the 
sample boundaries. The decoupled ly = ±n cases with 
conserved are adiabatically connected, once is not 
conserved, to the ordinary insulator for even n and to the 
topological insulator for odd n. A review of how these 
two cases emerge as the only possibilities in 2D follows 
in Section H. 

It is not obvious at first glance how to generalize the 
topological insulator phase to finite, noncrystalline sys- 
tems, rather than band structures, as when the parame- 
ters of the Hamiltonian H = Hq + Hso + Hfj are drawn 
from a random distribution. The first approach was in 
terms of a spin Chern number^ similar to the Chern in- 
teger in finite IQHE systems, but there is now agreement 
that for a clean band structure the only in variants are 
of Z2 type, rather than integer type.l^^l^JJiSI Two equiva- 
lent definitions of the appropriate Z2 invariant for a finite 
disordered system, in the simple case when the disorder 
splits all degeneracies other than Kramers degeneracies, 
are as follows (the full definition is given and compared 
to previous work in the following section) . The finite sys- 
tem can be considered as a unit "supercell" of a large 2D 
lattice. A large, finite supercell gives many bands, but 
each pair of bands connected by time reversal (Kramers 
pair) can be assigned its own Z2 invariant The phase 
of the supercell system, if the Fermi level lies in a gap, is 
then identified by adding up all the invariants (mod 2). 
Alternately, a direct definition of the phase in the finite 
system can be given that is related to the notion of "Z2 
pumping" Real charge is pumped as the fiux through 
the system is taken from to hc/2e (half the usual flux 
quantum that appears in IQHE pumping) ; we show that 
in the topological insulator, any pumping cycle, prop- 
erly defined, pumps an odd number of electron charges, 
while for the ordinary insulator any cycle pumps an even 
number of charges. 

We implement this definition numerically using an ex- 
plicit algorithm introduced by Fukui and HatsugaJ^ for 
computing Z2 topological invariants on a Brillouin zone. 
The topological insulator phase is robust to disorder: 
while different realizations of disorder assign different 
"Chern parities" to individual subbands, it is found that 
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the total for occupied subbands is always "odd" for a 
wide range of parameters, which in our definition indi- 
cates a topological insulator. In the IQHE, a pair of 
bands of opposite Chern number can annihilate as the 
strength of disorder is increased; in the QSHE, two band 
pairs that both have odd Chern parity can annihilate, 
i.e., become two even-parity band pairs. If the topologi- 
cal insulator can be destroyed by band annihilation, then 
there are extended (i.e., topologically nontrivial) states 
with an arbitrarily small gap; it may be the case that 
for some range of parameters, there are extended states 
at the Fermi level even in the thermodynamic limit, in- 
dicating a metallic phase. In the IQHE, there is only a 
single energy with extended states rather than a range 
of energies, and hence no metallic phase. We find the 
phase diagram of the graphene model with on-site dis- 
order, and in the presence of non-zero Rashba coupling 
find evidence for a metallic phase intervening between 
ordinary and topological insulators. 

Recent work by Obuse et al.^^ obtains a phase di- 
agram and critical exponents using a network model 
for the spin quantum Hall effect that is similar to the 
Chalker-Coddington network modeP^ for the IQHE (see 
also Onoda et alP^ for a quasi- ID study of localization 
in the Kane-Mele Hamiltonian with disorder). Our re- 
sults on the phase diagram are consistent with these 
works, although our method is unable to generate large 
enough system sizes to confirm the exponents found for 
the phase transitions. To understand how the network 
and Chern-parity approaches complement each other, 
consider the integer quantum Hall effect (IQHE): while 
the phenomenological network approach to the IQHE is 
valuable both to find the critical indices precisely and 
to identify the minimal necessary elements of a theory 
for the transition, Chern-number studies remain impor- 
tant f or studie s of effects such as the floating of extended 
states,^ ^^ * ^^ ! where knowledge of the topological prop- 
erties of a state is required. The network model gives 
more accurate information about the phase transitions 
but, if only the localization length is probed, does not 
distinguish the different phases in bulk. A more techni- 
cal difference between the two approaches is discussed at 
the end of Section II. 

Section II reviews how the topological insulator phase 
in perfect crystals arises from a parity-valued topological 
invariant of the band structure, similar to the TKNN^° 
integers in the IQHE. It then gives two mathematically 
equivalent definitions of the topological insulator phase 
in disordered systems based on Chern parity. One def- 
inition simply considers a finite disordered system as a 
supercell of an infinite lattice system, while the other 
is based upon closed charge pumping cycles driven by 
application of fiux to a finite periodic system. Section 
III reviews the Fukui-Hatsugai algorithrrP^ adapted to 
numerical computation of these invariants in disordered 
systems, then computes the phase diagram of the Kane- 
Mele graphene model^-* with on-site disorder. The con- 
clusions of our study for general 2D disordered systems 



are summarized in a short Section IV. Whil e the re is a 
three-dimensional version of the QSHE^-'^ 2II23] ^^iq,^ jg 
less directly connected to the IQHE and has interesting 
localization behavior, we will restrict our attention to 2D 
except for some comments in the final section. 

II. CHERN PARITIES FOR DISORDERED 
NONINTERACTING ELECTRON SYSTEMS 

A. The definition of the Z2 invariant in clean 
systems 

We review one definition of the Z2 invariant of a band 
pair in a 2D band structure, then explain its generaliza- 
tion to noncrystalline systems. With the Hamiltonian 

H = Ho + Hso + Hr (5) 

defined in equations ([2]), ([s]), and Q, or any periodic, 
single-electron Hamiltonian, a Berry connection can 
be defined on the Brillouin zone (BZ) from the periodic 
part u{k) of a Bloch state ipu — M(k)e''' ''. For a single 
nondegenerate band, we define Jij{k) ~ — Vfc|uj) for 
band j. This Berry connection can be used to understand 
how a commensurate magnetic field in a two-dimensional 
lattice system, which allows definition of Bloch states on 
the (magnetic) Brillouin zone, leads to the IQHE. The 
TKNN integer^^ in the IQHE can be written as integrals 
over the BZ of the Berry field strength = (Vfc x A.)z, 

rij = — — / Tjd^k 
271' Jbz 



du 


du \ 


1 du 
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\dky 
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The TKNN integers are integer-valued topological in- 
variants: the TKNN integer of a band cannot change as 
long as the band remains nondegenerate. When bands 
touch, only the total TKNN integer of the bands is well- 
defined ,'231 and the IQHE state of a gapped system is given 
by the sum of TKNN integers for occupied bands. How- 
ever, all these integers vanish in a time-reversal-invariant 
band structure (see, e.g., Ref . [TTt . Instead there is a Z2 
invariant associated with a band pair in a time-reversal- 
invariant 2D Fermi system.!^ Time-reversal invariance re- 
quires that 

en{-k)e-^ ^nik), (?) 

where 7i(fc) is the Bloch Hamiltonian and O = —icTyK 
is the action of time reversal {K performs complex con- 
jugation and (jy, the usual Pauli matrix, acts on spin 
indices). 

With time-reversal breaking, as in the IQHE, it is rea- 
sonable to ignore band degeneracies, but time-reversal 
invariance forces the single electron energies to be doubly 
degenerate at certain time-reversal-invariant momenta 
(see, for example, Ref. . These degeneracies are known 
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as "Kramers degeneracies". The single-band connection 
A.j will not be globally defined. So long as each Kramers 
pair remains separated by gaps in the energy spectrum 
from all others, the appropriate generalization is 



A;j = -i{{Uji\Vk\Ujl) + (Uj2|Vfc|Uj2)) 

= —i tr u^^kUj. 



(8) 



In the second equation, Uj — {uji,Uj2) where the 
Bloch functions are viewed as column vectors, i.e., Uj 
is a matrix. This compact notation follows Fukui and 
HatsugaP^ and makes it simple to include more than 
one band pair simply by adding more columns to u. Fi- 
nally, we again define a field strength associated with the 
potential, T — (Vfc x A)z- 

From these quantities, Fu and Kand^ give the follow- 
ing formula for the Z2 topological invariant in terms of 
the Bloch functions of the clean system: 



D 



1 

2n 



dk A 



d{EBZ) 



EBZ 



mod 2. (9) 



The notation EBZ stands for Effective Brillouin Zone, 21 
which describes one half of the Brillouin zone together 
with appropriate boundary conditions. Since the BZ is 
a torus, the EBZ can be viewed as a cylinder, and its 
boundary 9(EBZ) as two circles, as in Fig.[2Fb). While T 
is gauge-invariant, ^is not, and different (time-reversal- 
invariant) gauges can change the sum of boundary inte- 
gral by an even amount. Heuristically, (|9| states that 
integrating the field strength of a Kramers pair of bands 
over half the Brillouin zone gives a topological invari- 
ant in the same way as integrating the field strength of 
a single band over the whole zone. However, a gauge- 
dependent boundary term must be added to make the 
integral an integer. As mentioned, this boundary term is 
ambiguous (via gauge changes) up to an even integer, so 
only the parity D mod 2 is well-defined. 
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(b) 
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A direct proof of the existence and Z2 nature of the 
topological invariant for multiple bands without intro- 
ducing gauge-dependent quantities can be obtaineci^] by 
considering maps, from the EBZ to the space of Bloch 
Hamiltonians, that are consistent with time-reversal, fol- 
lowing work on the IQHE by Avron et al.'^* The Berry 
field strength can be written in terms of the (gauge- 
invariant) projection operator onto the band pair rather 
than the (gauge-dependent) wave functions. In this ap- 
proach, the ambiguity by an even integer that is crucial 
to obtain a Z2 rather than Z invariant corresponds to 
the many different ways in which the circles that form 
EBZ boundaries in Fig. ^h) can be contracted to make 
the EBZ into a sphere. The boundary integrals in ^ 
just calculate the contribution to the Chern number from 
these "contractions." On this sphere, Chern integers are 
well-defined for each nondegenerate band pair, but the 
different ways of contracting the boundaries cause the 
resulting integers to differ by even numbers. An explicit 
numerical implementation of the Fu-Kane formula ^ 
was given by Fukui and Hatsugal^ and will be reviewed 
in Section III. 

The invariant D classifies the topology of band pairs; 
each pair carries a value of D = or 1. The model of 
graphcnc we study has four bands (two sites per unit 
cell, two spin states per site), so two band pairs Di and 
D2, a nd th e total band structure must have D — 
mod 2^^^ so that Di = D2- Hence there are two 
phases, characterized by the value oi D = Di, say. Now, 
if at some value of parameters the upper and lower band 
pairs meet somewhere in the BZ, D is not well-defined, 
so there can be a phase transition. The two pairs with 
D = I annihilate into two D = pairs when the bands 
collide as the parameters are varied. Note that, just as 
total Chern number is conserved when bands collide in 
the IQHE, total Z2 is conserved here. 

At half filling, these phases are both insulating, since 
there is a gap between the second and third bands; the 
semimetal that appears when the gap closes marks the 
phase boundary. What happens when disorder is added? 
We now explain how the definition of the topological in- 
sulator generalizes to disordered systems. Just as TKNN 
integers of a band structure give rise to Chern integers 
in a disordered system of finite size, the Z2 invariants of 
band pairs become "Chern parities" . 



The topological insulator phase for Slater 
determinants via Chern parity 



FIG. 2: (a) A two-dimensional Brillouin zone; note that 
any such Brillouin zone, including that for graphene, can be 
smoothly deformed to a torus. The labeled points are time- 
reversal-invariant momenta, (b) The effective Brillouin zone 
(EBZ). The horizontal lines on the boundary circles 9(EBZ) 
connect time-reversal-conjugate points, where the Hamiltoni- 
ans are related by time reversal and so cannot be specified 
independently. 



The TKNN integers generalize in the presence of dis- 
order and interactions to the Chern number of the many- 
particle wavefunction.^® A natural question is how disor- 
der and interactions modify the Z2 invariants of band 
pairs in spin-orbit-coupled 2D band structures. All 
derivations of the Z2 invariants of clean systems depend 
on Fermi statistics in some way: for example, the exis- 
tence of Kramers degeneracies and the related fact that 
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the time-reversal operator squares to — 1 both depend on 
Fermi statistics. A many-fermion wavefunction describ- 
ing an even number of fermions does not behave in the 
same way as single-fermion wavefunctions under time- 
reversal. Hence, given only the many-fermion wavefunc- 
tion, it does not seem likely that there is a generalization 
of the Z2 invariant. 

However, for the particular case of many-fermion wave- 
functions that are single Slater determinants of single- 
particle wavefunctions, the invariant can be generalized 
as we now show. While the assumption of a single Slater 
determinant limits the treatment of interactions to the 
Hartree-Fock level, there is no requirement that the wave- 
functions in the Slater determinant be Bloch states. As 
a result, the topological insulator and QSHE can be de- 
fined and studied for any disorder strength. 

Niu et al?^ and Avron and Seileif^S showed that for dis- 
ordered quantum Hall systems there exists a generaliza- 
tion of the TKNN invariant defined for clean systems.-- 
They introduce generalized periodic boundary condi- 
tions and find an invariant Chern number, similar in 
form to the TKNN invariant, on the space of bound- 
ary phases. Consider a finite system of noninteracting 
electrons with boundary conditions that are periodic up 
to phases 0a;, ^y, as shown in Fig. [l] this is equivalent 
to putting magnetic fluxes '^x,y = 4'x,y^a/'^''^ through 
the two noncontractible circles on the torus (<i>o — hc/e 
is the magnetic flux quantum). As motivation, think of 
the flnite system as a (possibly very large) unit cell of a 
lattice system. Then in order to determine the phase of 
this lattice system, instead of integrating over k to do the 
integrals in the Fu-Kane formula (|9| , our strategy will be 
to integrate over the boundary phases, which introduce 
offsets to the wave vectors. We first carry out this proce- 
dure, show that it reproduces the band-structure result 
for clean systems, and then discuss a physical picture and 
its relation to previous definitions. 

Consider the single-particle wavefunctions of a lattice 
Hamiltonian such as the graphene model on a finite lat- 
tice of size Lx X Ly (see Fig. fTl). Instead of the physical 
boundary conditions ip{x + Ll. y) — i^{x) for a single- 
particle wavefunction V', introduce the boundary phases, 
or "twists", 4> = {(f>x,4'y) via 

+ L.^) = e'"^-V'(a;), i^{x + Ly) = e^'^^^x). (10) 

In the special case of a clean system, this shifts each al- 
lowed value of wave vector fc by A = {(f)x/Lx,(f)y/Ly). 
Because it is simpler numerically to work with single- 
valued wavefunctions than the multiple-valued ones of 
equation (10 1, and because the standard formahsm of 



the Berry connection assumes the parameters to be ex- 
plicit in the form of the Hamiltonian rather than living in 
the boundary conditions (i.e., the definition of the Hilbert 
space) , we want to transfer the twist angles to the Hamil- 
tonian. This is done by making a unitary transformation 
to 



Note that for the special case (j)x.y = 27r the complex 
phase in (11) is single- valued as we translate the coor- 



dinates around the torus, so that this is just a gauge 
transformation; if we take the phases in the square 
[— TT, vr] X [— tTjTt], opposite boundaries are identified un- 
der smooth gauge transformations and the twist space is 
a torus. Under the change of basis (111, Kane and Mele's 
model Hamiltonian^ H = Hq + Hso + Hn (see Section I) 
becomes (suppressing spinor indices) 



H H{(f>) = ^ c/ [t + iXr{s X 4), 



(12) 



where again A = {(px/ Lx,4>y/ Ly), dij is still the vector 
i j, and we have added a random term for on-site 
disorder, Wi, drawn from the Gaussian distribution of 
zero mean and standard deviation a^. It is now clear 
that under time reversal, 



eH{-(t})e-'^ ^n{(t)), 



(13) 



since the extra phase factors in Ti. change sign under com- 
plex conjugation. Since this directly parallels equation 
(0, the Z2 invariant D of the Brillouin zone passes di- 
rectly to twist space: 



1 

2tt 



d4> A 



d(ETZ) 



ETZ 



mod 2, 



(14) 

with ETZ for Effective Twist Zone, i.e., ETZ = {</> | < 
(j^x < TT, — TT < 0J, < tt}. Note that there is an indepen- 
dent Chern parity for each Kramers-degenerate band 
pair separated from the rest of the spectrum by a gap at 
all 4>: for such an isolated pair, A. and T are defined as 
in ([8| and after, with k ^ <p and u ^ x- 

In order to make contact with the band-structure defi- 
nition, we note that if there is no disorder in the Hamilto- 
nian (i.e., (Jtu = 0), there are discrete translational sym- 
metries within the x Ly supercell that induce addi- 
tional non-Kramers degeneracies at some points in twist 
space. With such degeneracies only the total Chern par- 
ity of all the degenerate states is well-defined. For ex- 
ample, in the clean graphene model, there are only two 
separated sets of states, even though the dimension of 
the Bloch Hamiltonian increases as L^ and Ly increase. 



X{x) 



ip{x). 



(11) 



t Note that this means "H is not generically time-reversal invariant. 
Indeed, there are only four boundary conditions that respect 
time-reversal, namely those for which i/j picks up a real phase 
upon translation around each cycle of the honeycomb lattice. 
These four correspond to the TRIM of the clean system in the 
calculation of D. 
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Each of these sets has its own Chern parity, consistent 
with the counting of Z2 invariants in the band struc- 
ture (note that the two invariants are not independent 
because of a zero sum rule) 1^ As in the band-structure 
case, the phase of a physical system is determined by the 
sum of Chern parities for occupied band pairs. Disorder, 
as discussed in the following section, breaks all degenera- 
cies resulting from translational invariance, leaving only 
separated band pairs, each of which has its own Chern 
parity. We now discuss to what extent Chern parities can 
be connected to observable quantities in a finite system. 



C. Charge pumping cycles in 
time-reversal-invariant systems 

The above definition was motivated by thinking of a 
finite system as a (possibly large) unit cell of an infinite 
periodic lattice. The finite system was defined to be in 
the topological insulator phase according to whether the 
infinite lattice is in the topological insulator phase. Since 
the original identification of this phase was by numeri- 
cal observation of stable edge states in finite systems, it 
seems clear that there should be a direct way to detect 
the topological insulator in the finite system. 

The total Chern number in a finite IQHE system 
can be interpreted as measuring the number of charges 
pumped when the flux through one noncontractible circle 
on the torus increases adiabatically by one flux quantum. 
Briefly, one of the boundary phases corresponds to this 
driving flux, and the average over the other can be shown 
to yield the pumped charge.l^ The idea of "Z2 pumping" 
suggested by Fu and Kan^^^ is the following: in a finite 
cylinder with boundaries, the operation of increasing the 
phase (j)x (in the periodic direction, around the cylinder) 
from to TT, corresponding to a magnetic flux of one-half 
flux quantum through the cylinder, has the following ef- 
fect in the topological insulator. The values (t)x = Q and 
(j)x = 1^ are special because, unlike general values, they 
are consistent with time-reversal invariance. At these 
special fluxes there are gapless states at the Fermi level 
that are localized near the edges because of the bulk gap. 
Fermi statistics requires that these states lie in Kramers 
doublets. If at zero flux, the Kramers doublet at one 
edge is partially occupied (has one state occupied), then 
the operation of changing the flux changes its occupancy 
to either double or zero occupancy. Since total charge 
is conserved, this requires a flow of Z2 from one bound- 
ary to the other. Note that in addition to the average 
over applied flux, there is effectively an average over the 
spatial extent of the edge, because the edge states are 
extended. 

We now give an alternate definition of the topological 
insulator in terms of cyclic pumping of ordinary charge. 
This definition is mathematically equivalent to the defini- 
tion of |IIB| based on treating the finite system as a super- 
cell. In order to describe a closed pumping cycle, we need 
to add a second stage to the Fu-Kane process of increas- 



ing boundary phase 4>x from to tt: although both these 
phases are consistent with time-reversal invariance, phys- 
ical properties, like the occupancy of an edge doublet, are 
not identical at these different values of (t>x (even at the 
same (j)y). The only requirement on the second stage is 
essentially that it return the system to its original state 
without applying a time-reversal-breaking flux. The sec- 
ond stage gives a closed pumping cycle that returns the 
system to its original Hamiltonian, which means that the 
amount of charge pumped in a specified cycle is a well- 
defined integer quantity. 

Although the number of charges pumped is dependent 
not only on the first stage but on the second stage as well, 
whether this is an even or odd number is entirely deter- 
mined by the first stage, as we now show. A closed pump- 
ing cycle is shown in Fig. [sj (As before, we discuss the 
cycle in terms of variable Hamiltonians in a fixed Hilbert 
space rather than in terms of changing the boundary con- 
ditions on Hilbert space.) The original physical system's 
ETZ is the first stage of the cycle: the Hamiltonians are 
functions of (t)x from to tt and 4>y from — tt to tt, with 
time-reversal constraints that act on the boundary circles 
ai (l)x — and (f>x — tt. If (t>x takes one of these values, 
then the Hamiltonian at (j)y is time-reversal conjugate to 
that at (j)y — —4>y The second stage can be any contin- 
uous change of the Hamiltonians that takes the (f>x = 
system back to the (j)x = system and always satisfies the 
conjugacy condition between <j)y and (p'y. This is the key 
difference between the second stage and the first stage: 
for intermediate values < (f>x < tt , there is no such con- 
jugacy condition. The physical interpretation is that the 
second stage should be possible without introducing flux 
through the first noncontractible circle. 

Now the torus shown in Fig.|3]has one Chern integer for 
each isolated band pair. Summing over occupied bands 
gives the amount of charge pumped in the cycle. Al- 
though this integer charge depends on the second stage, 
its parity is solely determined by the first stage, i.e., the 
physical system. In particular, for the ordinary insula- 
tor there is some closed cycle that pumps zero charge, 
while for the topological insulator there is some closed 
cycle that pumps unit charge. These results follow from 
the same proof as for the band structure case in Ref. [TT] 
one shows that the differences in resulting Chern integers 
between any two second stages are even. The pumping 
definition gives some physical intuition for the "contrac- 
tions" introduced there; instead of contracting the EBZ 
to a sphere, here the ETZ is contracted to a torus by 
adding the second stage. The technical reason that these 
two constructions are equivalent is that, since the ap- 
propriate spaces of Hamiltonians are contractible (i.e., 
have TTi = 0), the two closed manifolds, the torus and 
the sphere, both have the same topological invariants, 
namely integer- valued Chern numbers. 

The topological insulator in disordered systems has 
been studied previously by locating the transition be- 
tween topological and ordinary insulators as a point or 
region where the localization length of single-electron 
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0x = O 
= B 



, = JT 

c c 



A A r 

stage I: Stage II: 

insert flux, breaking T complete the cycle 

FIG. 3: Graphical representation of charge pumping cycle for 
Ghern parities. The first stage takes place on the ETZ (as in 
equation ( |14[ )), and the flux (fix increases adiabatically from 
to TT. In the second stage the Hamiltonian at {<j)x ~ it, <j)y) is 
adiabatically transported through the space of Hamiltonians 
to return to the Hamiltonian at {(j>x = 0, The difference 
between the second stage and the first is that at every step 
of the second stage, the Hamiltonians obey the time-reversal 
conditions required at = or (px = tt. The bold lines in- 
dicate paths along which all Hamiltonians are time-reversal 
invariant, and the disk with horizontal lines indicates, as be- 
fore, how pairs of points in the second stage are related by 
time-reversal. 



eigenstates diverges. The existence of the topological in- 
sulator is inferred from the existence of this transition 
region (or alternately from the edge states in the topolog- 
ical insulator phase). In principle this approach is differ- 
ent from ours, in that our definition probes the existence 
not just of extended states but specifically of extended 
states that contribute to the pumping of charge, or alter- 
nately can be expected to give rise to edge states. The 
same distinction arises in the quantum Hall effect, where 
looking for extended states of nonzero Chern number is a 
more direct probe of quantum Hall physics than consid- 
ering the inverse participation ratio, for example, which 
would detect extended states of zero Chern number in 
addition to topological states. However, as in the quan- 
tum Hall case, we find that the phase boundaries from 
our Chern-parity definition are consistent with thos e ob- 
tained from calculations of the localization length.'i^ESI 



III. GRAPHENE MODEL AND NUMERICS 



A. The phase diagram of the disordered graphene 
model 



conduct in the thermodynamic limit. In the presence of 
a magnetic field, as in the IQHE, there are isolated en- 
ergies with extended states, but no finite-width range of 
energies with extended states, and hence no true metal- 
lic phase. Hikami et al. showecP^ that disordered two- 
dimensional systems with spin-orbit coupling can never- 
theless support a metallic phase, referred to as a "sym- 
plectic metal". In our case, any metallic phase would 
presumably appear in a region around the parameter set 
that closes the (clean) gap, as depicted schematically in 
Fig.g 

We also know that at A/{ = the z component of spin 
is a good quantum number (s^ commutes with H), so 
the system reduces to two copies of the Haldane model,'^ 
which has a quantum Hall plateau transition with no 
metallic phase. 



B. Lattice Implementation 

For numerical work we use the algorithm of Fukui and 
Hatsugai, which we review here.^"^ The formula (|9| for D 
requires a gauge choice for the Hamiltonian eigenstates 
at each on the two boundaries of the half-torus < 
4>x < TT, — TT < (j)y < TT. That is, the "field strength" T is 
gauge invariant, but the gauge potential is not. The 
eigenstates form Kramers pairs related by time reversal, 
and the gauge choice must respect this constraint. 

Now, at the time-reversal invariant points cf> = (0,0), 
(0, tt), (tt, 0), and (tt, tt), the solid points in Fig. p\ the 
spectrum is degenerate, with two states at each energy. 
The gauge condition requires that Q interchange the two 
with a phase factor e*"/^ (so that 0^ = — 1 as required 
for single-fermion states) . Numerical diagonalization will 
not, in general, return eigenvectors that obey this condi- 
tion, but we can force them to do so as follows: choosing 
one of the two members of each Kramers pair at energy 
£2n-i = £2n and calling that vector x^n-i, we discard 
the other and replace it by 



X2„ 



Bx2n-1- 



(15) 



On the rest of the boundary, eigenvectors can be chosen 



freely on < 
takes 



< vr. On —7r<(j)y<0 the algorithm 



Xn(-0) = 0X«(0)- 



(16) 



In summary, the algorithm leaves alone the results of 
numerical diagonalization at all the gray points in Fig. [5] 
and by hand enforces the gauge condition on the rest of 
d(ETZ).i 



It is useful to review some general expectations be- 
fore applying the definitions of the previous section to 
study a disordered version of the graphene model. In two- 
dimensional systems with time-reversal invariance and no 
spin-orbit coupling, even very weak disorder will localize 
electron wavefunctions, so that these systems do not ever 



t The points {(f>x,TT) and (</>2,,— tt) are not physically distinct; as 
pointed out earlier, a gauge transformation relates ^{(px, —n) 
to ?i(</>a;, tt) as given, for example. We therefore impose 
'H{if>x, —tt) = 'H(<fix, t) in the calculation. Alternatively, we could 
introduce boundary terms along (jiy = tt, — tt to compensate for 
the discrepancy. 
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FIG. 4: Schematic piiase diagram, (a) Phase diagram for a clean system, with fixed t and Ago 7^ (after Kane and Mel^. 
(b) Again a clean system, now with fixed t and A„ 7^ 0. (c) The expected form of the phase diagram at nonzero disorder (we 
run all simulations at fixed A„). The phase boundary in (b) opens up into a metallic phase, closing only when A_r = 0, where 
there should be an IQHE transition. 
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FIG. 5: Twist Space. The bold lines indicate the bound- 
aries of the "Effective Twist Zone", the region we integrate 
(or sum) over to calculate the Chern parity. The arrows in- 
dicate the direction to perform the sum over the boundary 
terms, and the lattice sites in gray indicate those for which the 
Hamiltonian eigenvectors are independently specified. That 
is, time reversal symmetry determines the eigenvectors on the 
white sites once those at the gray sites are found. 



With the eigenstates fixed at each point on the ETZ, 
we follow Fukui and Hatsugai and construct C/(l) parallel 
transporters on the links, as 



9x 
\9x 



detxH(f')x{4^ + x) (17) 



and Uy similarly. Like m in ([8|, x is a matrix built from 
occupied state vectors, and x translates by one link in 
the (px direction. In the continuum limit g should ap- 
proach a pure phase, but for non-zero lattice spacing it 
will in general have \g\ < 1, since the occupied subspace 
of interest will not embed in the total Hilbert space in 
the same way at every lattice point. 



In the end, the only retained information will be the 
variation of the relative phases (hence the definition of 
U), which can be captured by choosing a lattice con- 
stant so small that the phase field varies slowly over one 
link. However, the scale of variation will presumably dif- 
fer for different disorder realizations, and we would like 
a way to diagnose this and throw out those realizations 
for which fast variation makes the calculation unreliable. 
The phase is periodic, and if it were to wind through 2t: 
over the distance of one link the algorithm would miss 
this fact, so the relative phase will not provide a good 
diagnostic. As a proxy, we choose to cull out disorder re- 
alizations that result in small determinants in most sim- 
ulations, since a small overlap between adjacent occupied 
eigenspaces indicates rapid variation. Of course, this fil- 
tering could introduce a selection bias into the results; 
these effects are within the statistical uncertainty of our 
analysis (in particular, within the error bars of Fig.[9|. 

Associated with the transporter on each link is a 
gauge potential A^^y = logUx^y. This A is pure imag- 
inary, and the logarithm is defined to return the branch 
A/i £ (— 7r,7r). Associated with the transport around 
each plaquette is a fiux 

^^(0) =logC/,(</,)C/y(0 + i)[/-i(0 + y)[/-i(0), (18) 



again satisfying F/i £ (— 7r,7r). With these definitions in 
hand, the lattice Z2 invariant corresponding to (14 1 
is 



Dr 



I 



\ediETZ) DeETZ 



mod 2. (19) 



Of course, the sum over the boundary should have the 
same orientation as the corresponding contour integral 
m ([W]); in Fig.|5] this means the sum on the left bound- 
ary should carry a minus sign, following the arrows. As 
mentioned previously, this formalism has the desirable 
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property of guaranteeing that Dl = or 1 , but using too 
coarse a mesh can return the wrong value. 



Numerical Results 



As noted after equation (14), the number of nondegen 



erate Kramers pairs in a disordered system will generi- 
cally be extensive, and we can use equation ( |19| to calcu- 
late the Chern parity of each pair separately. Fig.[6]shows 
the results of such calculations on a 4 x 6 lattice to present 
a picture of the three phases: normal insulator, symplec- 
tic metal, and topological insulator. In all phases, there 
are Kramers pairs with ZJ^ = 1 in the lower half of the 
energy spectrum, as indicated by the bars. However, in 
the normal insulator (Aso well below the transition, or 
gap closure) there are an even number of such pairs in 
the occupied half of all realizations, so that the overall 
Chern parity is even. The presence of a very small num- 
ber of realizations (2 of 215 for a 6 x 8 lattice) with odd 
parity indicates either the tail of the disorder-broadened 
transition or, more likely, that the weak filter applied 
for these simulations failed to catch all realizations with 
rapidly varying link variables. 
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FIG. 6: Distribution of Chern parities for band pairs. The bar 
heights represent the fraction of disorder realizations (out of 
~ 200 trials) that have a given number # of band pairs with 
Dl = 1 in the occupied (half-filled) subspace. Those with an 
even number # will have an overall Dl = 0, those with an 
odd number will have overall Dl = 1- All these simulations 
were done with t = — 1, Afl = A„ = l, and = 0.3. Read- 
ing across. Ago increases and the system transitions from all 
realizations having even parity at small Ago to odd parity at 
large Xso- Reading down, doubling the system size doubles 
the total number of Kramers pairs and roughly doubles the 
number of Z2-odd pairs. 

When Xso is large, almost all realizations have an odd 
number of Z2-odd Kramers pairs (as with the small-Aso 
case, 2 realizations do not follow the rule), and in the 
region near the transition of the clean system there are 



instances of both types. The presence of disorder causes 
the "gap" to close at different values of Xso for different 
realizations, and also at different energies. The latter 
fact means that extended states are present throughout 
a finite spread of energies, while the former means that 
the metallic state exists over a finite region of parameter 
space. 

For an integer quantum Hall system, Yang and Bhatifi^ 
have shown how to extract the localization length expo- 
nent V from such calculations, in their case the Chern 
numbers of Landau sublevels. Specifically, sublevels with 
non-zero Chern number contain extended states, which 
should only occur at isolated energies in the IQHE. 
Therefore, the number Nc of such sublevels should de- 
crease as the system size Ns increases, and in fact {Nc) (x 
n} ^I"^^ ^ where () indicates an average over disorder re- 
alizations. A similar approach for the QSH system here 
should reveal (A^'d) cx A^s, where Nd is the number of 
Z2-odd Kramers pairs, since we expect a stable metallic 
band of energies in the thermodynamic limit (as observed 
by Obuse et a/.'^and Onoda et al^. With enough data 
and large enough systems, finite-size-induced broadening 
of the edges of this band should also make v accessi- 
ble via a subleading term in the scaling. We find that 
larger system sizes require a finer mesh in twist space for 
these Kramers-pair-resolved simulations to return stable 
results, so that the requirements quickly outstrip our re- 
sources. Nevertheless, comparison of the two rows in 
Fig. [6] indicates that the location of the mean roughly 
doubles. This is what would be expected for the mid- 
dle case given the above considerations; the mean for the 
topological insulator (the right-hand panels in Fig. [6]) 
should grow more slowly, as in the case studied by Yang 
and Bhatt. 

The total phase of the system, given by the total Chern 
parity, is more relevant to possible measurements than 
the Chern parity of each Kramers pair. The former 
maintains its meaning if we consider the ground state 
wave function of the many-electron system, formed as a 
Slater determinant of the single-electron states we use 
here. There is also a computational benefit to calculat- 
ing Z?0 for the whole occupied subspace rather than for 
individual pairs, as well — at larger system sizes, and 
also at stronger disorder, the link variables for the half- 
filled subspace vary much more slowly than those for the 
individual Kramers pairs, making the calculation more 
robust. For these reasons the remaining plots in this pa- 
per depict only the Chern parity of the half-filled system. 

To show that a metallic region of non-zero extent in pa- 
rameter space exists in the thermodynamic limit, we need 
to verify that the mixed-phase region does not shrink to 
zero as we increase the system size. Figure |7][a) shows 
that as we make the system larger, the transition region 
certainly does not get narrower, and in fact the largest 
system size seems to have the broadest transition. 

We can quantify the scaling of the metallic region's 
width with system size by assuming a simple one- 
parameter scaling form for the curves in Fig. [7];a) and 
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true). In that case, the variance of the distribution can 
be estimated as cr^ cx p(l — where p is the fraction 
of disorder reahzations returning Z?^ = 1 , and then l/cr^ 
is an appropriate weight for the statistical test. The er- 
ror bars in Figs. [7] and |8] are also assigned based on a 
binomial model of the data.l^ 

By contrast, at = the Hamiltonian ([5) reduces to 
two copies of the Haldane model and so should exhibit the 
quantum Hall plateau transition, which looks like a step 
function at zero temperature in the thermodynamic limit. 
In Fig. |8ja) the width of the transition region shrinks as 
the system size grows, consistent with the prediction. 
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FIG. 7: (a) At finite A_r, the "metallic" region persists as 
the system size grows, and even broadens in the case shown 
here {t — —1, Xr = = ~ We identify the metallic 
region as those values of Xso for which some, but not all, 
disorder realizations have Dl = 1. As explained in the text, 
the fit to the simulation data (shown only for 6x8 systems 
here) has the form (tanha + l)/2, with a = m{Xso — ^*so)- 
The error bars represent 95% confidence intervals assuming 
a binomial distribution of outcomes for each Xso- (b) The 
scaling collapse of the data in (a), based on the best fit m 
and A Jo) for each system size. 
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defining the width of the curve to be proportional to the 
reciprocal of the maximum slope: widths 1/slope. With 
sufficient data, one could expand the approximately lin- 
ear region near the middle of the transition in a power 
series with a few coefficients as fit parameters. Since our 
simulation data are limited, we opt instead to assume 
the form (tanh a{s) + 1) /2, which has roughly the right 
shape. If a = m(Xgo ~ -^so) ' f ti^n m is exactly the maxi- 
mum slope we want. Figure|7jb) plots the data versus the 
best-fit scaling variable a for each system size; the points 
appear to fill out a smooth curve, justifying the scaling 
hypothesis. The best fit for m and Acn is determined 
by minimizing a weighted statistic.'^ In particular, 
we assume that the results (Dl = ±1) of independent 
simulations at a fixed parameter set are distributed bi- 
nomially, and that for each system size there are the same 
number of trials at each value of Xso (which is roughly 



FIG. 8: (a) At Xr = 0, the metallic region gets narrower as 
the system size increases. The fit is to a 6 x 6 system and 
Xn = 0; as in Fig. [t] t — —1, A„ = cr„ = 1. (b) The scaling 
variable a again comes from fitting to a tanh, and error bars 
represent 95% confidence intervals. 

More quantitatively, PruiskeiP^l has shown in a renor- 
malization group framework that the functional form of 
the crossover for the IQHE looks like 

piL,B)^f{a), a (X L^^^iB - B*), (20) 

where L is the linear size of the finite sample, B is the ap- 
plied magnetic field, and v is again the localization length 
exponent. The function p could be either the longitudi- 
nal or transverse conductivity in the IQHE. Therefore, 
p{oo,B* ± e) = /(±oo), i.e., the transition is sharp in 
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the thermodynamic limit. In our system the analogous 
parameter to B is Xso- in the Haldane model, the spin- 
orbit coupling breaks time-reversal invariance locally, like 
B does in the IQHE. The width of the transition region 
in B is governed by the way the Landau level energies 
respond to changing B, and the width of the transition 
region in our model (for fixed Ay) is determined by the re- 
sponse of the gap to changing Xso ■ Since both responses 
are linear, we expect that the appropriate scaling variable 
will be a cx L^/^iXso - A^o)- 

This form would allow us to extract the exponent i/ 
from the scaling of the maximum slope with system size 
for large systems (there are corrections at small system 
sizes). Again making the fit to a tanh described above. 
Fig. [sjb) shows the scaling form, and Fig. [9] shows the 
scaling of width with linear system size. In particular, 
a regression gives l/iy = 0.78 ± 0.03, to be compared 
with the accepted value of 1/j^ w 0.42. That is far from 
good agreement, but the observed scaling should not be 
taken as implying a new universality class. (For refer- 
ence, the network-model work by Obuse et al^ found 
1/u ^ 0.37, and Onoda et aZ.li^ recently found a value 
l/i^ Ki 0.63) First, there should be finite-size corrections 
to the simple scaling assumed for the small systems con- 
sidered here. Second, the scaling form is in principle dif- 
ferent for different geometries, and the simulations were 
done for systems of varying aspect ratio (1 to 1.5). Never- 
theless, it is clear that the qualitative behavior at A;? = 
is as expected, showing no metallic phase, and the behav- 
ior at A/j = 1 is consistent with the presence of a metallic 
region. 
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FIG. 9: Width of the metallic transition region as a function 
of linear size L. The width is determined by the reciprocal of 
the maximum slope of the simulation data in Figures [7] and 
[8j and the linear size is taken to be the square root of the 
number of sites. 



Finally, by varying A r and noting the Xso values that 
mark the edges of the transition region for each Xn, 
we can map out the phase diagram of the Hamiltonian 
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FIG. 10: Approximate width of metallic region for a 4 x 6 
lattice (fixed X^ = = The dashed curves indicate the 
parameter values at which 98% and 2% of the disorder real- 
izations have Dl ~ 1- This underestimates the true width 
of the "metallic" region but hopefully avoids some amount 
of the inevitable error due to small system size. This dia- 
gram should offer a reasonable approximation to the thermo- 
dynamic (infinite-system) phase diagram away from A_r = 0. 
At Xr — 0, there should be no metallic phase, but a sharp 
transition between the two insulating phases in the thermo- 
dynamic limit. Given the plots in Fig. [7] it appears that finite 
size effects also reduce the width of the metallic phase at large 

Xr. 



([5]), as in Fig. 10 The widths obtained this way are in 
rough agreement with the scaling analysis outlined above, 
which returns the behavior of the width and not its nor- 
malization. As the results in Figures [7] and [8] show, this 
phase diagram will overestimate the width of the metallic 
phase at A^ = 0, which is really zero. Together, these 
simulations confirm the expectation of Fig. |4] within the 
accuracy of our computational methods. 



IV. SUMMARY 



Previous worlP ^° l ^^ l '^^ l defined a Z2 topological in- 
variant in infinite lattices that is similar to the TKNN 
invariant for the integer quantum Hall effect.!^ In disor- 
dered systems with boundaries, Fu and Kan^^ defined 
a topological invariant in terms of pumping of the oc- 
cupancy of Kramers-degenerate edge states. We have 
given a definition of a topological invariant valid for dis- 
ordered systems without boundary, i.e., without appeal 
to edge states. The "Chern parity" can be thought of as 
describing either a finite system with boundary phases 
or an arbitrarily large supercell in an infinite lattice sys- 
tem with well defined wavevector. A physical effect of 
Chern parity is that it determines whether the amount 
of charge pumped in a certain type of closed pumping 
cycle is even or odd. The "ordinary" and "topological" 
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insulator phases can be distinguished by this invariant as 
long as many-body effects do not prevent description of 
the ground state as a Slater determinant. Chern parity in 
the spin quantum Hall effect is the natural generalization 
of Chern number in the integer quantum Hall effect. 

In a disordered system, the only degeneracies expected 
to survive are Kramers degeneracies at time-reversal in- 
variant values of the boundary phases. In this case, each 
pair of states related by Kramers degeneracies can be as- 
signed its own Chern parity, and the overall Chern par- 
ity of all occupied state pairs determines the observable 
phase. The lattice algorithm for Z2 topological invari- 
ants laid out by Fukui and HatsugaP^ allows numerical 
identification of the topological insulator phase in disor- 
dered QSH systems. Implementing this algorithm for the 
specific graphene model Hamiltonian of Kane and Mel^ 
with added on-site disorder, we observe the ordinary and 
topological insulator phases in simulations. While the 
number of "odd" pairs (state pairs with odd Chern par- 
ity) varies with the disorder realization, there is an even 
number of odd pairs in the ordinary insulator, and an 
odd number of odd pairs in the topological insulator. 

We find that a metallic phase opens up between the two 



insulating phases for generic spin-orbit coupling. This 
agrees with the prediction of Hikami et al^ that spin- 
orbit coupling can protect a 2D metallic phase from disor- 
der, and confirms the simulation results of Onoda et al^ 
and Obuse et al}"^ The methods in this paper could in 
principle be used to study the three-dimensional case and 
confirm the argument in Ref. that only one of the four 
invariants of a band structure^^ ^is stable to disorder. 
While there is now strong evidenc^^^ that the phase tran- 
sitions in the 2D QSHE are, except for special points, in 
the previously studied symplectic metal-insulator class, 
there is as yet no numerical study of the phase transi- 
tions in three-dimensional topological insulators. 
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